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ABSTRACT 

We have performed A'^-body simulations of star clusters orbiting in a spherically symmetric smooth 
galactic potential. The model clusters cover a range of initial half-mass radii and orbital eccentricities 
in order to test the historical assumption that the tidal radius of a cluster is imposed at perigalacticon. 
The traditional assumption for globular clusters is that since the internal relaxation time is larger than 
its orbital period, the cluster is tidally stripped at perigalacticon. Instead, our simulations show that 
a cluster with an eccentric orbit does not need to fully relax in order to expand. After a perigalactic 
pass, a cluster re-captures previously unbound stars, and the tidal shock at perigalacticon has the 
effect of energizing inner region stars to larger orbits. Therefore, instead of the limiting radius being 
imposed at perigalacticon, it more nearly traces the instantaneous tidal radius of the cluster at any 
point in the orbit. We present a numerical correction factor to theoretical tidal radii calculated at 
perigalacticon which takes into consideration both the orbital eccentricity and current orbital phase 
of the cluster. 

Subject headings: galaxies: kinematics and dynamics - globular clusters: general 



1. INTRODUCTION 

Theoretical calculations of the radius of a star cluster, 
or tidal radius, usually assume that the gravitational field 
of the host galaxy regulates cluster size. In most previous 
treatments, for simplicity it is further assumed that the 
gravitational field in which the cluster orbits is constant, 
i.e. the cluster has a circular o rbit in a spherically sym- 
metr ic galactic potential (e.e. Ivon Hoerneriil957l IK ma 
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2003: iBinnev fc Tremaind 120081: iBertin fc Va rri 200? 
The tidal radius is then assumed to be equal to the Jacobi 
radius (rj), the distance beyond which the acceleration 
a star feels towards the galaxy center is greater than the 
acceleration it feels towards the cluster center, and the 
star is able to esc ape. First-order tid al theory determines 
the tidal radius (Ivon Hoerneil [19571) via: 
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(1) 



where Rgc is the galactocentric distance of the cluster, 
Mci is cluster mass, and Mg is the mass of the galaxy 
(assumed in early studies to be a point mass). 

For a cluster with a non-circular orbit, the fact that 
the tidal field is no longer static makes a n ana lytic ap- 
proach very difficult (see iRenaud et al.l ()2011h for an- 
other approach). Historically it has been assumed that 
for a globular cluster on an eccentric orbit, its tidal ra- 
dius is imposed at perigalacticon {Rp) where the tidal 
field of the host galaxy is the strongest. This ass ump- 
tion was in i tially suggested by Ivon Hoernerl (|1957( ) and 
later iKinel (jl962l ). and follows from the fact that the 
internal relaxation time (ir^) of the cluster is greater 
than its orbital period for almost all observed globular 
clusters. Therefore after stars outside the tidal radius 
at perigalacticon escape, the cluster would not be able 
to relax and expand before it returns to perigalacticon. 
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Thus in Equation |l] Rgc is usually replaced with Rp to 
calculat e the tidal radius of a cluster with an eccentric or- 
bit fe.g.'Innanen. Harris, fc Webbink''1983':'Fall fc ZhangI 
2001; Read ct al. 2006; Webb, Sills, fc Harris 2012). 

However, recent studies are showing with increasingly 
strong evidence that the actual sizes of observed clus- 
ters are not imposed at perigalacticon in this simple 
way. The actual size of an observed cluster is known 
as its limiting radius r^,, which ma rks the point where 
the c luster density falls to zero (|Binnev fc Tremaind 
|2008[) . Using the solved orbits of 15 Galactic glob- 
ular clusters, lOdenkirchen et all (|1997t ) demonstrated 
that cluster limiting radii are not solely depe ndent on 
periga lactic distance. iBrosche. Odenkirchen. fc Geffcrt] 
|l99y) suggested some sort of orbit- averaged tidal ra- 
dius is more appropriate when predicting limiting radii. 
Even with the orbits of an additi onal 29 Milky Way 
globular clust ers currently known (|Dinescu et al.l IT999I : 
Casetti-Dinesc u~et al.l l2007t ) , there is still no clear re- 
lationship between limiting r adii and perigalact i c dis- 
tance, and the conclusions of lOdenkirchen et al.l (|1997l ) 
stiU hold. 

With the Galactic potential as given bv I Johnston et ahl 
(1995) (the same Galactic potential Casetti-Dinescu et 
al. used to solve cluster orbits), we calculated the theo- 
retical tidal radius rt at perigalacticon (that is, the Ja- 
cobi radius at i? = i?p) of each of the 44 Galactic clusters 
with known orbit s. The rt value s were calculated with 
the formalism of iBertin fc Varri! (|2008l ) (as outlined in 
Section 3.0). The main uncertainty in the theoretical 
tidal radius is due t o the uncertainty in R^ quoted i n 
iDinescu et all (|1999D and ICasetti-Dinescu et all (|2007D . 
To compare theory and observations for individual clus- 
ters, we take clus ter limiting radii as determined from 
direct iKinsJ (|1966l ) model fits r^ to the observed cluster 
profiles as listed in|llarris (1996) (2010 edition). We then 
calculate the ratio of the difference (r^ — rt) between ob- 
served and theoretically predicted values to their average 
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Fig. 1. — Ratio of difference between observed (rj.) and theoret- 
ical (rt) tidal radius at perigalacticon to the average {(rt + rj,)/2) 
versus perigalactic distance for Galactic globular clusters. 

{{^k+i't)/'2). An uncertainty of 10% was assigned to val- 
ues of r/c. If theory and observations are in agreement, 
the ratio will be approximately zero. Clusters which have 
a ratio greater than zero will be ones which overfill their 
predicted tidal radius, while clusters with ratios less than 
zero will be tidally under-filling. The comparison be- 
tween theory and observations is illustrated in Figure [1] 

Whil e it is not expect ed that all clusters are tidally 
filling (jGieles et al.ll2010 f). the fact that the majority of 
clusters appear to be tidally overfilling is a strong sig- 
nal that the simple assumptions built into the model 
need investigation. The known pr esence of tidal tails 
around obse rved (e.g. [O dcnkirchcn c t al.ll200lD and sim- 
ulated (e.g . Mont uori et al.i i2007i: iKiipper et all 120121: 
iLane et al . 2012) globular clusters is not sufficient to ex- 
plain the cases of apparent overfilling. The observed 
values are determined from King-model profile fits that 
are heavily dominated by the inner populations of stars, 
out to a few half-light radii. In almost all cases the 
extremely low densities of stars in the tails, at or be- 
yond the nominal tidal radiu s , exer t little leverage on 
these fits. Furthermore, iKingI (|1966l ) models are known 
to underestimate cluster sizes in general as they require 
a sharp tidal cutoff whic h is not observed in all clusters 
(|McLaughlin & van der Marcl 2005 ). Only i n the most 
extreme cases (e.g. Pal 5 (Oden kirchen et al.ll2003| )). will 
large and extended tidal tails influence model fits to the 
cluster surface brightness profile. 

Theoretical calculations may also underestimate be- 
cause the assumption that the tidal radius is imposed at 
perigalacticon implies that the shape of the tidal field 
and the cluster orbit are not important. Hence factors 
such as tidal heating and disk shocking due to a varying 
symmetric galactic potential are not taken into consid- 
eration. However, despite these potential inaccuracies 
we still expect to see some sort of correlation between 
Tfc and perigalactic distance, if tidal radii are imposed 
at perigalacticon. The Milky Way cluster data therefore 



suggest that something is wrong with basic tidal theory . 

Recent A^-body si mulat i ons by IKiipper et al.l (|2010l ) , 
find that their fitted iKind ()1962t) radius was better rep- 
resented by the time averaged mean tidal radius of the 
cluster and not the perigalactic tid al radius. In a later 
study on the structure of tidal tails. iKiipper et an()2012l ) 
found that while stars outside the tidal radius as cal- 
culated at perigalacticon will likely become unbound at 
perigalacticon, some are able to be re-captured by the 
cluster as it moves away from perigalacticon and the in- 
stantaneous tidal radius of the cluster increases. That 
is, the limiting radius of a cluster will be greater than 
the tidal radius calculated at perigalacticon. This dis- 
crepancy is expected to be amplified for clusters on very 
eccentric orbits, as they make quick perigalactic passes 
and spend the majority of their time ne ar apogalacti- 
con (R g). In fact, iV-body simulations bv lMadrid et al.l 
(|2012D suggest that the half-mass radius of a globular 
cluster is more likely imposed at Ra than Rp. 

The purpose of this study is to explore more thoroughly 
the influence of orbital eccentricity on cluster size. Model 
A^-body clusters with different initial half-mass radii are 
evolved from zero age to 10 Gyr over a range of orbital 
eccentricities in the disk of a Milky Way-like potential. 
The models and their initial conditions are described in 
Section 2. In Section 3 we focus on the influence of or- 
bital eccentricity on the mass (Af ), half-mass radius {r„i), 
limiting radius, and tidal radius of each model cluster 
over time. In Section 4 we explore the influence of ini- 
tial cluster half-mass radius on our results. In Section 
5 we discuss the results of all our A^-body models and 
the infiuence of orbital eccentricity not only on cluster 
radii, but on individual stars within the cluster as well. 
Based on our findings, in Section 6 we suggest a correc- 
tion factor that can be applied to the perigalactic tidal 
radius of a cluster to better match its observed limit- 
ing radius. This correction factor is then applied to the 
Galactic globular clusters shown in Figure [T] Finally in 
Section 7 we summarize our conclusions. 

2. THE MODELS 

We use the NB0DY6 direct A^-body code (|Aarseth I 
l2003f) to evolve model star clusters from zero age to a 
Hubble time, over a range of both initial cluster half- 
mass radii and orbital eccentricity. All models in this 
study begin with 48000 single stars and 2000 binaries. 

As long as we use one particle to represent one real 
star, clusters of this size correspond physically to either 
a very massive open cluster, or a low-mass globular clus- 
ter in the Milky Way. Ideally, we would like to follow 
more "average" globular clusters, which will typically 
have 200,000 stars or more. However, to this point in the 
history of the subject, direct N-body integrations with 
N-values that large have only been car ried out for spe- 
cial, specific purposes (for example, see iHurlev fc Sharal 
(20_12) for a 200,000-star simulation in the tidal field 
of a point-mass Galacti c potential and a circular or- 
bit; IZonoozi et all ()2011[ ) for N < 100,000-star simu- 
lations d irected at modell i ng Pa l 14, again on a circular 
orbit; or iHeggie fc Gierszj ()2009[ ) for a 105,000-star sim- 
ulation of NGC 6397 over 1 Gyr). Encouragingly, how- 
ever, these high-N models generally confirm the trends 
obtained from earlier, small-N simulations (see Hurley et 
al. for additional discussion). 
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The purpose of our suite of 50,000-star models is in- 
stead to survey a wide parameter space of initial cluster 
half-mass radii and orbital type. Our chosen ranges (see 
below) match real Milky Way star clusters moving in 
a realistic time-varying potential. Eventually, with ad- 
vances in computational capabilities this type of survey 
work can be extended to higher N values (up to 10^) 
that will cover almost the entire known range of globu- 
lar clusters. However, as will be seen below, the models 
summarized here already prove to be highly informative 
in revealing important physical effects of orbital eccen- 
tricity on the internal dynamical evolution in a direct 
way that does not rely on analytical approximations. 

Since we are only concerned with the influence of or- 
bital eccentricity on clusters of different initial half mass 
radii, our choice of initial parameters such as cluster 
metallicity, the stellar initial mass function (IMF), and 
binary fraction are of little consequence as long as they 
remain consistent between models. However, we note 
that binary fractions of a few per cen t are typical for 
globular clusters (e.g. lDavis et al.|[2008[ ). 

The masses of single stars are drawn from a 
iKroupa. Tout, fc GilmoTi (|1993D IMF between 0.1 and 
30 Mq. For binary stars, the masses of two randomly se- 
lected single stars are combined to equal the total mass 
of the binary, with the primary and secondary masses 
determined by a mass-ratio randomly drawn from a uni- 
form distribution. The initial total mass of each model 
is 3 X 10'* Mq. The initi al period of each binary is d rawn 
from the distribution of lDuauennov fc Mavoii (|1991[ ) and 
their orbital ecce ntricities are a ssumed to follow a ther- 
mal distribution (jHeggie II1975I) . All stars were given a 
metallicity oi Z = 0.001. The initial positions and ve- 
locities of the stars are generated based on a Plum mer 
density profile ([Plummer l ligillAarseth et al.l[l97l . We 
note that the Plummer model extends to an infinite ra- 
dius so we impose a cut-off at 10 to avoid rare cases 
of large cluster-centric distances. A description of the al- 
gorithms for st ellar and binary evolution can be found in 
IHurlevI (12008^ . 

The model clusters orbit in a three-dimensional Galac- 
tic potential, which cons ists of a point-mass bulge, a 
iMivamoto fc Nagail (|1975f ) disk (with a = 4.5 kpc and 
h = 0.5 kpc), and a logarithmic halo potential. The 
combined mass profiles of all three potentials result in 
a circular velocity of 220 km/s at a galactocentric dis- 
tance of 8.5 kpc. The bulge and dis k have masses o f 
1.5x10*° and 5 x 10 |°M(7) respectively (iXue et al.ll2008l) 
lAarseth I (|2003[ ) and iPraagman. Hurley, fc Poweil (|2oTo[ ) 
describe the incorporation of the Galactic potential into 
NB0DY6. All models were set to orbit in the plane of the 
disk such that a cluster on a circular orbit experiences a 
static tidal field, and will not be subject to factors such 
as tidal heating or disk shocking. 

For the purposes of our study, the only parameters 
which are important and change from model to model are 
initial cluster half-mass radius, initial cluster position, 
and initial cluster velocity which determine the shape of 
the orbit. Our first models were for clusters with orbital 
eccentricities of 0.25, 0.5, 0.75, and 0.9, where eccentric- 



TABLE 1 
Model Input Parameters 



ity is defined as e = 



All of these have the same 



Ra+Rp 

perigalactic distances of 6 kpc and initial half-mass radii 
rm,i of 6 pc. These clusters are located at perigalac- 
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ticon at time zero. For comparison purposes a model 
was simulated with a circular orbit at perigalacticon (6 
kpc) and four more with circular orbits at the apogalactic 
distance of each eccentric cluster (10 kpc, 18 kpc, 43 kpc, 
and 104 kpc). These simulations allow us to directly com- 
pare the properties of a globular cluster on an eccentric 
orbit to clusters on circular orbits at both perigalacticon 
and apogalacticon. 

Our set of models also included re-simulations of the 
cluster with an orbital eccentricity of 0.5, and the corre- 
sponding e = perigalactic and apogalactic simulations, 
but with initial half-mass radii of 4pc, 2pc, 1 pc, and 
0.5 pc. This range allows us to study the influence of 
initial cluster half-mass radius. Both sets of models are 
summarized in Table [TJ Model names are based on or- 
bital eccentricity (e.g. e05), circular radius at apogalac- 
ticon (e.g. rl8), and initial half mass radius (rm6). 
Hence a model cluster with a perigalactic distance of 
6 kpc, apogalactic distance of 18 kpc (orbital eccentric- 
ity of 0.5), and an initial half-mass radius of 6pc would 
be labeled e05rl8rm6. 

3. INFLUENCE OF ORBITAL ECCENTRICITY 

The first portion of this study will focus solely on mod- 
els with an initial equal to 6 pc, with the only dif- 
ference between each model being their Galactic orbits. 
However, we first need to determine whether any given 
star is bound to the cluster. In a cluster-centric coordi- 
nate system, we define the x-axis as pointing away from 
the galactic center, the y-axis pointing in the direction 
of motion of the cluster, and the z-axis pointing perpen- 
dicular to the orbital plane. In this coordinate system 
the energy of an individual star can be written as: 
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Fig. 2. — Snapshot of a model cluster on a circular orbit at 6 
kpc after 32 orbits (5680 Myr). Bound stars are marked in red. 



where the second term is the potential energy due to the 
remaining N-1 stars in the simulation, each with mass 
rrii and located a distance from the star. The third 
term is the tidal potential with Vl^ equal to the orbital 
frequency of the cluster. Here u is a dimensionless pos- 
itive coefficient defin ed below in Equation [51 The tidal 
potential, taken from lBertin fc Varril ()2008[ ). results in a 
stretching of the cluster in the x-direction, no change in 
the y-direction, and a compression in the z-direction. If 
the resultant energy is less than zero the star is bound, 
otherwise it is considered unbound. 

We considered additional criteria for determining 
whether a star is bound or unbound in addition to 
the energy calculation. Other studies have invoked a 
distance cutoff such that the stars' cluster-centric dis- 
tance must be greater than the cluster perigalactic or 
instantaneous tidal radius for it to be unbound (e.g. 
iTakahashi &: Baumgardtl [20121 ). It has also been sug- 
gested that a star's velocity plays a role in whether or 
not it can be considered unbound (e.g. iKiipper et alJ 
[2OT0l [2O12). However, we found these additional criteria 
did not change any of the results found in Section 3 as 
they only effected a small percentage of simulated stars. 
Therefore, we only require that a star's energy as given 
by Equation [2] to be greater than zero for the star to be 
considered unbound. 

Figure [2] shows a model cluster at a representative 
timestep. The tidal tails formed by escaping stars are 
clearly visible in our simulations. The densely populated 
spherical collection of stars marked in red are those that 
satisfy our boundedness criterion and are considered clus- 
ter members. The unbound stars that appear close to the 
centre of the cluster are foreground stars with a large z 
coordinate and are simply projected onto the cluster in 
the x-y plane. These tails have no effect on our determi- 
nation of theoretical tidal radii or observed limiting radii 
as we only consider stars that are bound to the cluster. 



3.1. Mass 

We first wish to study how orbital eccentricity influ- 
ences the total mass, or more specifically the mass loss, 
of a star cluster over time, which then plays an important 
role in determining tidal radii (r^ oc Af a ). The total mass 
of stars bound to the cluster in each model is illustrated 
in Figure |3l In general, mass loss is due to stellar evolu- 
tion, evaporation due to two-body interactions, and tidal 
stripping. For clusters with circular orbits, in all cases 
the apogalactic cluster loses mass at a lower rate than 
the perigalactic case as a result of less tidal stripping. 

The e=0.25 case loses mass at almost the same rate as 
if it had a circular orbit at its perigalactic distance, but 
the final mass is still notably larger than the perigalactic 
case. Since eccentric clusters spend the majority of their 
time away from perigalacticon, they too will be subject 
to less tidal stripping than an ideal cluster that spends 
all its time at perigalacticon. As eccentricity increases, 
the mass-loss profile shifts further away from the peri- 
galactic case and closer to a cluster with a circular orbit 
at apogalacticon. 

At higher eccentricities the mass no longer smoothly 
decreases, in contrast to the circular orbit cases. Instead 
periodic fluctuations are present. The minima of these 
fluctuations correspond to perigalactic passes, where the 
rapid increase in tidal fleld strength results in episodes of 
significant mass loss. These fiuctuations suggest that a 
greater change in tidal field strength between apogalac- 
ticon and perigalacticon results in stars gaining more en- 
ergy at or near perigalacticon. 

Especially interesting in the lower right panel of Fig- 
ure [3] is the fact that once a cluster undergoes significant 
mass loss during a perigalactic pass, the cluster starts 
to regain mass before resuming its mean mass loss rate. 
In these intervals just after perigalacticon, the cluster is 
re-capturing some of the stars which were previously un- 
bound. Stated differently, many of the stars that were 
formally unbound at Rp drift away slowly enough that 
they are recaptured as the cluster moves back outward 
and its instantaneous tidal radius expands again. Fur- 
thermore, these fluctuations suggest that while a peri- 
galactic pass does have a strong effect on an eccentric 
cluster, the cluster cannot be treated as if it had a circu- 
lar orbit at R„. These re s ults a re in agreement with the 
findings of IKiipper et al.l ()2012D discussed earlier. 

3.2. Half-mass Radius 

We next consider how orbital eccentricity can influence 
the half-mass radius Vm of bound stars within a globular 
cluster. It should be noted that r,„ is not the same as 
the half-light radius (also known as the effective ra- 
dius), which is a directly observable parameter. In our 
simulations, the half-mass radius is always slightly larger 
than the half-light radius. 

The results of our simulations are illustrated in Figure 
m The initial increase in during the flrst ^ 2000 Myr 
in all cases is driven by two-body relaxation and stellar 
evolution mass-loss. However once the cluster is relaxed, 
tidal stripping becomes the dominant dynamical process. 

The Tm profiles of the apogalactic cases in the lower 
panels do not begin to decrease after 2000 Myr, but in- 
stead continue to increase up to 10 Gyr. As discussed in 
the next section this trend is due to the fact that these 
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Fig. 3. — The evolution of total cluster mass over time. The red 
lines correspond to models with orbital eccentricities as labelled 
in each panel. In each plot, the lower black line corresponds to 
a cluster with a circular orbit at perigalacticon (6kpc), while the 
upper black line corresponds to a cluster with a circular orbit at 
apogalacticon. 



clusters are barely tidally filling, so can still expand and 
not be subject to tidal stripping. 

Similar to the results of Section 3.1, for low eccen- 
tricities the r„i profile of the eccentric model cluster is 
comparable to the circular orbit case at perigalacticon. 
Increasing eccentricity brings the profile closer to the 
apogalactic case on the average, but increasing eccen- 
tricity again results in sharp fluctuations in the pro- 
file which correspond to perigalactic passes. The trends 
shown in Figure |4] reveal what is perhaps the most strik- 
ing difference between clusters on static circular orbits 
(the classically assumed case) and ones on more realistic 
eccentric orbits. If cluster limiting radii are imposed at 
perigalacticon, we would expect the minima of the eccen- 
tric profile to be equal to the profile of a cluster 
orbiting at Rp. While a high-e cluster may briefly ex- 
pand after a perigalactic pass, the next perigalactic pass 
would restore to a size equal to the perigalactic case. 
But not even at perigalacticon does the of the eccen- 
tric cluster equal the perigalactic case. Instead, the time 
averaged is reflective of a tidal fleld w eaker than the 
field a t perigalacticon, in agreement with iKiipper et al.l 
(|20T0l) . 

After a perigalactic pass, Figure|4]illustrates again that 
the cluster is able to increase in size. Especially ap- 
parent in the lower right panel of Figure U is the fact 
that the cluster is able to increase to a size greater 
than its mean r^- Inspection of our iV-body models 
shows that this increase in size is due to a combina- 
tion of re-capturing s ome of the previously unbound stars 
(jKiipper et al.l[2013 ) and the stars in the inner region of 
the cluster gaining enough energy to move outward and 
repopulate the halo of the cluster. These statements are 
discussed in further detail in Section 5.0. 
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Fig. 4. — The evolution of half mass radius over time. The red 
lines correspond to models with orbital eccentricities as labelled 
in each panel. In each plot, the lower black line corresponds to 
a cluster with a circular orbit at perigalacticon (6kpc), while the 
upper black line corresponds to a cluster with a circular orbit at 
apogalacticon. 

3.3. Tidal and Limiting Radii 

The Jacobi radius represents a theoretical surface 
around a cluster past which a star cannot pass and still 
remain bound. The Jacobi radius allows for the calcu- 
lation of the instantaneous tidal radius of each model 
cluster as a function of time. To calculate the instan- 
taneous tidal radius of each model cluster, we require a 
derivation of cluster tidal radius which takes into con- 
sideration the tidal fleld of the host galaxy. Assuming 
only that the tidal fleld must be spheric ally symmetric, 
the theoretical tidal radius as derived by lBertin fc Varo 
((200a is: 



where f2, n and v are defined as: 



n'' = id^GiR)/dR)RjRp 



{dHG{R)/dR^) 



(3) 

(4) 
(5) 
(6) 



$G is the galactic potential, Rp is the perigalactic dis- 
tance, O is the orbital frequency of the cluster, k is the 
epicyclic frequency of the cluster at Rp, and u is a posi- 
tive dimensionless coefficient. Using the tidal field of the 
Milky Way discussed in Section 2.0 and the mass and 
galactocentric distance of the model clusters at each time 
step, we calculate the instantaneous tidal radius of each 
model. The results of these calculations are shown in Fig- 
ure [S] Here the instantaneous increases and decreases 
periodically along the orbit, but never quite reaches the 
perigalactic and apogalactic cases due to differences in 
mass loss rates among all three cases. 
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Fig. 5. — The evolution of the instantaneous tidal radius over 
time. The red lines correspond to models with orbital eccentricities 
as labelled in each panel. In each plot, the lower black line corre- 
sponds to a cluster with a circular orbit at perigalacticon (6kpc), 
while the upper black line corresponds to a cluster with a circular 
orbit at apogalacticon. 

Next we compare the tidal radius to the hmiting ra- 
dius. For a simulated cluster, since we know which stars 
are bound or unbound, we could call the limiting ra- 
dius of the cluster the distance to the farthest bound 
star, but this approach introduces some significant prob- 
lems. First, since cluster tidal radii are calculated for 
stars with circular prograde orbits, any star with a ret- 
rograde and/or eccentric orbit within its clu ster can re- 
main , bound beyond the nominal tidal radius (jRead et al.l 
I2006D . Second, any star in the process of escaping the 
cluster can reach large clustercentric distances before be- 
coming energetically unbound. Third, the stars along the 
y-axis of the cluster (the direction of motion) are unaf- 
fected by the tidal potential in Equation^ allowing them 
to also remain bound at larger clustercentric distances. 
These three issues cause the true limiting radius of the 
cluster to change dramatically from time-step to time- 
step. To gain a more stable indication of cluster size, 
we instead focus on the x-axis of the cluster, the axis 
along which the tidal radius is calculated, and define the 
limiting radius as the average x-coordinate of all stars 
with ||x|| > Tf. This calculation typically involves less 
than 1% of the total cluster population. While this is 
not the true limiting radius of the cluster and will al- 
ways be slightly larger than the true tidal radius, it acts 
as a tracer of the outer region that is less affected by 
individual extreme outliers. If a cluster is tidally over- 
filling, the limiting radius will still be significantly larger 
than the tidal radius. For a cluster that is tidally under- 
filling, the limiting radius is simply the distance to the 
outermost bound star. 

In Figure [6] we show this empirically determined r^ for 
each model as a function of time. For circular orbits, 
on average the limiting radius of the cluster decreases 
smoothly as a result of mass loss. For eccentric orbits. 
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Fig. 6. — The evolution of the limiting radius over time. The red 
lines correspond to models with orbital eccentricities as labelled 
in each panel. In each plot, the lower black line corresponds to 
a cluster with a circular orbit at perigalacticon {6kpc), while the 
upper black line corresponds to a cluster with a circular orbit at 
apogalacticon. 

the small fluctuations with perigalactic passes in Figures 
[3] and 2] are much more prevalent in Figure [51 Comparing 
Figure [5] to Figure [SJ the fluctuations in both figures 
indicate that the limiting radius behaves the same as the 
instantaneous tidal radius. 

It should be noted that the apogalactic cases for e — 
0.75 and 0.9 in Figure [6] are not smooth due to the fact 
that these clusters are barely tidally filling and their lim- 
iting radii are easily influenced by individual escaping 
stars. 

Directly comparing limiting radii in Figure |6] and tidal 
radii in Figure [5l all circular orbits have a relatively con- 
stant ratio at approximately ^ = 1.1. Since we expect 
r^ to slightly overestimate rt, a ratio of 1.1 suggests that 
these clusters are approximately tidally fllling. For ec- 
centric clusters the ratio is in general also 1.1, suggesting 
the clusters come close to fllling their instantaneous tidal 
radius at all times. Fluctuations in the ratio for the e = 
0.75 and 0.9 cases indicate that after a perigalactic pass 
the cluster is slightly tidally under-fllling and works to 
fill its instantaneous tidal radius on the way to apogalac- 
ticon. When travelling back in from apogalacticon to 
perigalacticon, the cluster will remain tidally filled and 
lose stars to tidal stripping as the instantaneous tidal 
radius shrinks. 

4. INFLUENCE OF INITIAL CLUSTER HALF-MASS 
RADIUS 

Up until this point we have only considered clusters 
with initial half-mass radii of 6 pc. This initial half-mass 
radius was chosen simply to ensure that the model clus- 
ters with Rp — 6kpc would be tidally filling. As seen 
in Figure S] this produces clusters with sizes at 10 Gyr 
ranging from 3 to 14 pc, which are larger than most (but 
not all) real globular clusters. In an attempt to produce 
Milky Way-like clusters which have a mean effective ra- 
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Fig. 7. — The evolution of half mass radius over time. The 
upper left, upper right, lower left, and lower right panels are for 
simulations with initial half mass radii of 4 pc, 2 pc, 1 pc, and 0.5 pc 
respectively. The red lines correspond to models with orbital eccen- 
tricities of 0.5, while the lower black lines correspond to a cluster 
with a circular orbit at perigalacticon (6 kpc) and the upper black 
lines correspond to a cluster with a circular orbit at apogalacticon 
(18 kpc). 

dius of 2.5 pc, we re-simulated the e0r6rni6, e05rl8rm6, 
and e0rl8rm6 models with initial half-mass radii of 4pc, 
2pc, 1 pc, and 0.5 pc. The results are shown in plots of 
half mass radius versus time in Figure [71 

The rm4 clusters closely resemble the rm6 clusters, 
with similar periodic fluctuations with perigalactic passes 
and the final values for the perigalactic, eccentric, and 
apogalactic cases. However, the rm4 clusters undergo 
smaller initial expansion due to two-body interactions 
than the rm6 clusters, and thus when outer region stars 
are being removed through tidal stripping, since the ma- 
jority of the mass is concentrated in the inner region, the 
mass-loss profile is less affected. 

This issue of initial size becomes even more significant 
in the rm2, rmf, and rm05 cases. The periodic fluc- 
tuations on the eccentric orbit are barely visible in the 
rm2 cluster, and non-existent in the rml and rm05 clus- 
ters. The rml cases are significantly smaller at 10 Gyr, 
approaching the typical ^ 2 - 3pc size that match the 
majority of real globular clusters. The rm05 models have 
completely dissolved by 10 Gyr. 

These small clusters are only tidally filling in the sense 
that two-body interactions have pushed some bound 
stars to orbits that take them out to the instantaneous 
tidal radius. With the majority of the bound stars lo- 
cated in the inner regions of the cluster, tidal stripping 
is not the dominant form of mass loss and the influence 
of the galactic potential and cluster orbit are minimized. 
Instead, stellar evolution and two-body interactions are 
the dominant forms of mass loss. These clusters would 
be cla ssified as "tidally unaffected" (ICarballo-Be llo et al.l 
[Ml ). For the rml case, the profiles of the perigalac- 
tic, eccentric, and apogalactic cases begin to split only 
after ^ 5 Gyr, when stars have been finally pushed to the 



outer regions of the cluster and tidal stripping is begin- 
ning to play an important role. While this is true, it is 
not due to two-body interactions but instead a result of 
core collapse. For the rm05 case, not even core collapse 
can push stars to the outer region of the cluster in order 
for tidal stripping to occur. In fact, the rm05 clusters all 
have the same r™ up until the complete evaporation of 
the cluster at approximately 7 Gyr. 

Producing Milky Way-like globular cluster effective 
radii of 2 - 3 pc for clusters on circular orbits at 6 kpc 
or greater appears to require initial size less than 
1 pc. The iV-body models reveal that either clusters 
originally are extremely compact and tidally unaffected, 
or present-day cluster orbits have changed significantly 
from the orbit along which the clusters originally formed. 
Some observational support for this view can be found 
in recen t measuremeiits of very young, massive clus- 
ters f e.g.lBastian et al."2008l. l2012t[iPortegies Zwart et al.l 
120101: [Marks fc K roupa 2010 ^ How ever, recent A^-body 
simulations by [Sippel et al.l (|2012 f) showed that and 
be very different because of stellar mass segrega- 
tion, and produce clusters with final values near 3 pc 
despite large rm.i- These issues will be explored in future 
studies. 

5. DISCUSSION 

A perigalactic pass has three effects on a globular clus- 
ter, which we illustrate in Figure [8] for the e=0.9 model 
e09rl04rm6. In this figure, we plot the energy per unit 
mass of individual stars (as per Equation[2l) as a function 
of radial distance from the cluster center at 9 points in 
the orbit. Beginning in Panel A of Figure [HI for a given 
time between apogalacticon and perigalacticon there are 
a few stars that are within close proximity of the clus- 
ter but remain unbound (marked in red). As the clus- 
ter moves towards Rp and the instantaneous tidal radius 
shrinks (Panel B), more and more stars become tem- 
porarily unbound. As predicted in Section 3.0, even stars 
in the inner region of the cluster are provided with a sig- 
nificant increase in energy by the tidal shock and can be- 
come unbound. Just after the cluster reaches perigalac- 
ticon (Panel C), a large number of stars are no longer 
bound to the cluster. As the cluster moves away from 
perigalacticon (Panels D to I): 

• some stars that became unbound escape the cluster 
(which causes the initial decrease in and r^); 

• some of the stars that are unbound in Panel C re- 
turn to energies below zero and are recaptured (see 
Figure [9]); 

• the tidal shock gives stars initially found in the in- 
ner region enough additional energy to move out- 
ward and fill the orbits vacated by stars which per- 
manently escaped the cluster (see Figure [TU)) . 

It is even possible for inner region stars to become tem- 
porarily or permanently unbound if they move outward 
at a rate faster than the instantaneous tidal radius in- 
creases. 

6. PREDICTING CLUSTER LIMITING RADII 

Now that we have shown that limiting radii are not 
imposed at perigalacticon, it is useful to know how to 
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Fig. 8. — The radial distance and energy of stars within model 
cluster e09rl04rm6 for different time steps. Beginning in Panels A 
and B the cluster is travelling towards perigalacticon. In Panel C 
the cluster has just left perigalacticon. In Panels D to I the cluster 
is moving away from perigalacticon. Bound stars are marked as 
black and unbound stars are marked as red. 




2000 4000 6000 8000 10* 2000 4000 6000 8000 10' 

'I'lme (Myr) Time (Myr) 

5x10' 




3x10' 




2000 4000 6000 8000 10* 

Time (Myr) 



2000 4000 6000 8000 10* 

Time (Myr) 



Fig. 9. — The total number of bound (red) and unbound (green) 
stars in a cluster as a function of time. The black line corresponds 
to the total number of stars and the orbital eccentricity of the 
model is labelled in each panel. For comparison purposes, the 
number of unbound stars has been increased by a factor of 10. 
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Fig. 10. — The evolution of the number of stars within the inner 
10 % Lagrangian radius over time. The red line correspond to 
the e=0.9 model (e09rl04b). The lower black line corresponds to 
a cluster with a circular orbit at perigalacticon (6kpc), while the 
upper black line corresponds to a cluster with a circular orbit at 
apogalacticon (104 kpc). The number of stars within the inner 10 
% Lagrangian radius will naturally decrease over time due to two- 
body encounters, however the eccentric case (red line) illustrates 
that with each perigalactic pass a significant number of stars move 
beyond the inner 10 % Lagrangian radius. 

calculate a meaningful number that predicts the limiting 
radius of a globular cluster on an eccentric orbit. As we 
saw in Figures [5] and [6l the limiting radius essentially 
traces the instantaneous tidal radius. However the ratio 
between cluster limiting radius and instantaneous tidal 
radius undergoes small periodic variations as a function 
of the location of a cluster along its orbit. 
Orbital phase is defined as: 



F : 



Rn 



R-n 



Ra Rp 



(7) 



such that the cluster has F = at perigalacticon and 
= 1 at apogalacticon. A cluster with a circular orbit 
will always have F = Q. The median limiting radius of 
the eccentric cluster as a function of orbital phase is then 
determined in order to calculate the ratio of the instan- 
taneous limiting radius to the limiting radius of the e — 
perigalactic case ( r[Xe-=o) ^ ' ^^^^ normalized by mass. 
This ratio is plotted as a function of phase in Figure [TTJ 
It is important to note that we have ignored the first 2000 
Myr of evolution for each model cluster when evapora- 
tion due to two-body relaxation is the dominant source 
of mass loss. 



For a given orbital eccentricity, 



''i-(e) 



changes almost 



ri,(e=0) 

linearly with phase F. It is interesting to note that we ob- 
served a second order effect that the rate at which the 
cluster expands is lower than the rate at which it con- 
tracts. When the cluster is moving away from perigalac- 
ticon, it works to fill its expanding tidal radius. Con- 
versely a cluster moving towards perigalacticon would 
be larger as it is always tidally filling on the way inward. 



9 




0.4 0.6 

Phase (F) 

Fig. 11. — The ratio of the mass normalized limiting radius of 
model clusters with eccentric orbits to the mass normalized lim- 
iting radius of a cluster with a circular orbit at perigalacticon as 
a function of orbital phase as defined in Equation [7] Error bars 
represent an uncertainty of Icr. 



TABLE 2 
Lines of Best Fit 



Orbital Eccentricity Slope Uncertainty 



0.25 
0.5 
0.75 
0.9 



0.512 
1.29 
3.37 
7.84 



0.007 
0.04 
0.07 
0.07 



The rate at which the hmiting radius of the cluster 
increases and decreases as a function of orbital phase 
is a strong function of orbital eccentricity. These rates 
were determined explicitly by finding the slopes of each 
line, where the y-intercept is forced to equal one. The 
slopes are listed in Table [21 along with the associated 
uncertainty (Icr). 

A smooth relationship between slope and eccentricity 
emerges, that can be fit with an exponential, and allows 
us to propose a purely analytical correction to the cal- 
culation of the tidal radius of a globular cluster. For 
a globular cluster on an orbit with eccentricity E, with 
a tidal radius at perigalacticon rt{Rp) and located at a 
phase F in its orbit, its limiting radius is equal to 



rL{F) =rt{Rp)il + a F e'' 



(8) 



where a = 0.17 ± 0.03 and 6 = 4.1 ± 0.2. Note that 
since this calculation involves a single cluster over the 
course of a single orbit, the mass normalization is no 
longer necessary as cluster mass will not have changed 
significantly over a fraction of one orbit. 

The next step is to simulate a larger suite of model clus- 
ters ranging in initial cluster mass and half-mass radii to 
determine if these parameters play a role in the correc- 
tion factor suggested above. However, regardless of the 
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Fig. 12. — Ratio of difference between fitted iKinel l l 19621 ) radius 
(rfe) and Hmiting radius (rj^) to the average of the two radii versus 
perigalactic distance for Galactic globular clusters. Limiting radii 
have been calculated based on the orbital eccentricity and phase 
of a cluster as given by Equation [8] 

influence of cluster mass or initial half-mass radius, all 
tidally affected simulations follow the rule that the lim- 
iting radius traces the instantaneous tidal radius rather 
well. Thus if full orbital information or phase F is un- 
known, the calculation of the instantaneous tidal radius 
is a reasonable estimate of the limiting radius of a cluster. 
For globular clusters in other galaxies, in which only their 
projected galactocentric distances are known, it may be 
possible to determine their theoretical tidal radius based 
on their present King radius rk ■ Future work will explore 
this possibility. 

6.1. Application to the Milky Way 

For many Milky Way globular clusters, their current 
galactocentric distance, orbital eccentricity and orbital 
phase are known. In Figure [12] the revised, fully cor- 
rected version of Figure [T] is illustrated, where rj;, is now 
the phase-corrected value from Equation [H We now 
see more tidally under-filling clusters and the scatter of 
points more nearly around zero. A stronger agreement 
between theory and observations emerges. Correcting 
for using a non-spherically symmetric potential in cal- 
culating tidal radii and improved methods for determin- 
ing observational limiting radii will likely strengthen this 
comparison further. 

7. CONCLUSIONS AND FUTURE WORK 

Globular clusters have been simulated with a range 
of both circular and eccentric orbits. After determining 
which stars are bound to the cluster at a given time, 
we show that while eccentric clusters undergo episodes 
of significant mass loss during a perigalactic pass, their 
time averaged mass loss rate reflects a tidal field less 
than the tidal field at perigalacticon. Additionally it 
was found that clusters are able to re-capture unbound 
stars after a perigalactic pass as their instantaneous tidal 
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radius increases. 

Second, we show that the half-mass radius of a globular 
cluster increases and decreases about a mean value over 
the course of an orbit. These fluctuations suggest that 
the perigalactic pass also has the effect of energizing inner 
region stars to larger orbits. Finally, we find that the 
limiting radius of a cluster traces its instantaneous tidal 
radius at all times. 

These findings argue against the historical assumption 
that globular cluster tidal radii, and by extension limit- 
ing radii, are imposed at perigalacticon for clusters that 
do not have circular orbits. While it remains true that 
the half-mass relaxation time is greater than one orbital 
period, the cluster does not need to fully relax in order to 
expand. The eccentric orbit introduces an effect of tidal 
shocking that is not experienced by clusters in a static 
potential (circular orbit). 

While the instantaneous tidal radius is a useful first ap- 
proximation of the limiting radius, we have proposed an 
analytically determined correction factor that is a func- 
tion of orbital eccentricity and phase. This correction 
leads to a much stronger agreement between th e pre- 
dicted limiting radii and observational iKin^ ()1962[ ) radii 
of Milky Way globular clusters. Future studies will ex- 
plore how the correction factor depends on initial mass or 



initial half-mass radius and how corrected limiting radii 
are related to King radii. 

Since the tidal field of the Milky Way is not spherically 
symmetric, correcting limiting radii based on eccentricity 
and orbital phase is not the final step. We still need to 
correct for orbital inclination to account for factors like 
disk shocking and tidal heating, which may reveal im- 
portant effects for the Milky Way and other disk galax- 
ies. However, the present results already have clear ap- 
plicability to elliptical galaxies, which have more nearly 
spherical potentials. We are currently investigating iV- 
body simulations in these directions. The ultimate goal 
is to be able to predict the limiting radius of any tidally 
affected globular cluster, given its orbit, galactocentric 
position and the galactic potential of the host galaxy. 
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